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1 Abstract 

Several methods of incorporating multi-dimensional ideas into algorithms for the solution of 
the Euler equations are presented and discussed. Three schemes are developed and tested; 
a scheme based on a downwind distribution, a scheme based on a rotated Riemann solver 
and a scheme based on a generalized Riemann solver. The schemes all show a marked 
improvement over first-order, grid-aligned upwind schemes, but the higher-order performance 
is less impressive. An outlook for the future of multi-dimensional upwind schemes is given. 


2 Introduction 


2.1 The ID Euler Equations 

The Euler equations of gasclynamics express the conservation of mass, momentum and energy 
for a continuous, non-conducting inviscid fluid. In one dimension, they may be written 
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where u is the ‘state vector,’ the vector of conserved quantities, 
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and f is the ‘flux vector,’ given by 
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The ideal gas relation 
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and the definition of total enthalpy 
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close the set of equations. 

The primitive variable form of the equations may also be written in vector form, giving 
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where u is the vector of primitive variables 
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and a is the Jacobian matrix 
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Travelling wave solutions to the one-dimensional Euler equations may be sought by sub- 
stituting 

u (x,t) = u (x — Xt) (3a) 

into Equation 2a, giving the eigenvalue problem 


(a — AI) <5u = 0, 


(3b) 


where I is the three by three identity matrix, and <5u is the amplitude of the traveling wave. 
The solution of this problem depends upon the eigenvalues and eigenvectors of a described 
below. 
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First Acoustic Wave 
The eigenvalue 

A — u — c (4a) 

and corresponding right eigenvector 


(4b) 


represent an acoustic disturbance. The wave travels with the fluid velocity rrunus the 
speed of sound. 





also represent an acoustic disturbance. The wave travels with the fluid velocity plus 
the speed of sound. 

The “characteristic variables” are given by 

dw = r _1 6u (5a) 
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Modern schemes for compressible flow are based on these travelling wave solutions to 
the one-dimensional Euler equations. The Euler equations in their conservative form (Equa- 
tion la) are integrated over cell number t, which extends from cell-face * — 1/2 to cell-face 
t + 1/2. This gives 


/■•'+J du . f'+l m , 

L m dxJr L Tx ix = 0 


+i df 
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Defining the cell-average state vector u and applying Gauss’ theorem gives 
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+ ( f *i " 
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where fj+i /2 = f(uL, u R ) is a “numerical flux function,” constructed so that the components 
of dw for which A is positive (corresponding to a wave travelling in the positive x direction) 
are backward differenced, and the components for which A is negative (corresponding to a 
wave moving in the negative x direction) are forward differenced. In the above, ul and ur 
are values of the state vector just to the left and right sides of the interface, respectively. 
These can be taken to be equal to the nearest cell center values, yielding a first-order scheme, 
or obtained by some interpolation of the cell-center values, yielding a higher-order scheme. 
An outline of several of the flux functions currently in use has been given by Van Leer et 
«/[!]• 

The schemes outlined in this paper are primarily based on the flux function of Roe [2]. 
For this flux function, the “Roe-average” of the states ul and Ur is defined, with 
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The interface flux is then defined as 
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with ft, A k and AV4 as given below: 
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Ai = u — c , (11a) 

A 2 = A, (Hb) 

A 3 = U + C. (11c) 

As with the quasilinear form, the r4’s can be interpreted as waves of strength A 14, propa- 
gating with speed A*,. Here, the subscript 1 denotes the first acoustic wave, the subscript 2 
denotes the entropy wave and the subscript 3 the second acoustic wave. 


2.2 The 2D Euler Equations 

In two dimensions, the Euler equations may be written as 

dU dF dG _ 
dt dx dy ’ 
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where U is the state vector, 
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and F and G are the flux vectors, given by 
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The ideal gas relation for the two-dimensional case is 

+ t , 2 j 
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The Euler equations may also be written in the quasilinear form 
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where A and B are the Jacobian matrices 
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and U is the vector of primitive variables 
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u = 


I 
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(13d) 


Travelling wave solutions to the two-dimensional Euler equations are of the form 

U(x,y,t) = \J (xk x + yK y - Xt) (14a) 

where k = K x e x + « y e y is a unit vector in the direction of propagation. Substituting this 
form of solution into Equation 13a gives the eigenvalue problem 

(A*. + B K y - Al) SU = 0 , (14b) 


where I is the four by four identity matrix, and is the amplitude of the traveling wave. 
This deceptively simple-looking equation is much more complex than its one-dimensional 
counterpart. 

The matrix R is made up of the right eigenvectors of the matrix Ak x + Bk v . The four 
eigenvectors and their corresponding eigenvalues are as follows. 


First Acoustic Wave 
The eigenvalue 

X = uk x + VK y — c (15a) 

and corresponding right eigenvector 


Ri = < 


—CK, a 
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(15b) 


represent an acoustic disturbance. The wave travels with the projected fluid velocity 
minus the speed of sound. 


Shear Wave 

The eigenvalue 


A = UK X + VKy 


(15c) 



and right eigenvector 
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represent a change in velocity direction, and can be thought of as a shear wave. The 
wave travels with the projected fluid velocity. 


Entropy Wave 

The same eigenvalue 


and its other right eigenvector 


A = UK X + VK y 



(!5e) 
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represent an entropy wave. The wave travels with the projection of the fluid velocity 
onto the direction vector k, as above. 


Second Acoustic Wave 
The eigenvalue 


A = UK X -f VKy + C 


(15*) 


and right eigenvector 

f , 1 
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represent an acoustic disturbance. The wave travels with the projected fluid velocity 
plus the speed of sound. 

The “characteristic variables” for a given value of ( k x ,k v ) may be computed, giving 

dw = R -1 £u (16a) 
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This, however, does not say anything about the choice of the direction vector k=(k x ,k v ). 
Time-dependent and steady solutions of the Euler equations may include waves travelling in 
particular directions. It is the way in which the direction vectors for these waves are chosen 
that distinguishes the various wave-like models for the Euler equations. 

Most modern schemes for the two-dimensional Euler equations are based on grid-contravariant 
directions. The conservative Euler equations (Equation 12a) are integrated over a cell 
with boundary dflij. This gives 


ff dU rr (OF dG\ 

IL, ~m dA + Jh, 


The second integral is converted to a line integral over the cell boundary: 

d\J 


ff ~zpdA + <f {Fdy — Gdx) = 0 . 
JJ n,,, at Jd fij,. 
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The first integral can be expressed in terms of the change in the average state U in the cell, 
and the line integral becomes a sum over the faces (say, four) of the cell: 


dU 


A- + Z(FAy-GAx) t = 0. 


t=i 


After introducing the cell-face length As, 

As 2 = Ax 2 + Ay ~ , 

Equation 17c can be written as 

^ + n >AS < 
at ( _ l 

where F n is the flux normal to the cell face 
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cosOg = Ay /As , (17g) 

sin 6 g = -Ax/ As, (17h) 

u± = u cosOg — v sin O g , (17i) 

U|| = u sin O g + v cos 6 g . (17j) 

Roe’s flux- difference splitting consists of writing the interface flux as 

F(U L , Ur) = \ (Fl + Fr) - ± £ |A t | AV*Rj, , (18) 


where the eigenvectors are 
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AV is given by 


and the wave speeds are 


AV = < 
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Aj = u x -c, 

A 2 = «_l , 

A 3 = u ± , 

A 4 = u x + c. 


For the two-dimensional case, 
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and p, v and h 0 are as before. 
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2.3 Wave Models for the 2D Euler Equations 

For a more general two-dimensional scheme, directions other than the grid-contravariant 
directions should be chosen. The direction vector k can not be determined directly from 
the equations of motion, however. In general, in multi-dimensional flow, information may 
propagate in an infinite number of directions. Thus a dominant direction, or directions, must 
be chosen for a discrete wave model, based on local flow values. 

There are basically two approaches to determining the propagation directions: 

1. an ad hoc choice of a dominant direction, or directions, based on physical considera- 
tions; 

2. using local flow gradients to “fit” a set of discrete waves to the residual, solving for the 
appropriate convection directions. 

The first approach is by far the simpler of the two. In it, dominant directions are chosen, and 
the residual is interpreted as a sum of waves moving in these directions. This approach is 
actually the one used in most current flow solvers, with the dominant directions taken to be 
the grid contravariant directions. This is clearly not physical. If, for instance, a stationary 
normal shock were to lie oblique to the grid, the difference between the pre- and post-shock 
states would be interpreted as a combination of a compression and a shear, instead of just 
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Figure 1: Representation of a Shock Oblique to the Grid 


as a compression (See Figure 1). If, however, the shock were studied in a frame normal to 
it, there would be no element of shear seen in the difference between the two states (See 
Figure 2). 

Clearly, then, the choice of dominant direction can have an effect on the solution. But 
what direction should be chosen? In the previous example, choosing the flow direction would 
have given the proper behavior, since the shock was normal to the flow direction. This would 
not be true for an oblique shock, however. For a non-curved oblique shock, directions that 
would give the proper behavior include: 

• the pressure gradient direction; 

• the density gradient direction; 

• the entropy gradient direction. 

Another direction of interest is the velocity difference direction, 

ur ~ ul 

ur — ui ) 2 + ( vr — vt ) 2 
vr-v l 

f\y — # ■ 7 ■ ■ "~ V " ■ = 

V («fl - ul ,) 2 + ( v/i — vl ) 2 

In this direction, the difference in the two states (Ur, Uj,) may be interpreted (based on 
the velocities) as either a, shock aligned normal to the velocity difference direction, or a shear 
aligned parallel to the velocity difference direction (See Figure 2). Not until the difference in 
pressure between the two states is taken into account can it be determined whether a shock 
(large pressure difference) or a shear (zero or small pressure change) exists. 

The second approach is more intricate, and has led thus far to two methods of calculating 
propagation directions. The first method is based on decomposition of local gradients, and 
is due to Roe [3]. The second method is based on an approximate diagonalization of the 
Euler equations, and is due to Hirsch et al [4] and Deconinck et al [5]. 
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(23a) 

(23b) 





Figure 2: Velocity Difference Direction 


In the first method, an a priori model for the number and type of waves is made, and 
then angles, strengths and speeds of the waves are determined from local data. The example 
carried out here will consisit of four acoustic waves, a shear wave and an entropy wave. The 
acoustic waves are taken to have unknown strengths ai, a 2 , <*3 and 04, and orientations 0, 
0 + t/2, 0 + 7 t, and 0 + 3 tt/ 2, respectively. The entropy wave is taken to have strength /3, and 
orientation <f>. The shear is represented by a uniform vorticity u> (See Figure 3). This model 
gives eight unknowns (a r , a 2 , <23, a 4 ,fl, u>, Q, <j>). The eight equations come from expressing 
the components of the local gradients of density, velocities and pressure in terms of these 
waves. The entropy wave, for example, contributes 

0R\ (cos <f>e x + sin <j>e y ) (24a) 

to the local gradient VU. The first acoustic wave contributes 

a x R 3 (cos $e x + sin 0e y ) (24b) 


where 


K x = COS 0 , 

Ky = sin 0 . 


The other acoustic-wave contributions are calculated similarly. The contribution of the 
vorticity is 

f 0 ) 
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Figure 3: Roe’s Wave Decomposition 


Summing the contributions of these waves gives the system of equations 


e*i cos 9 + a 2 cos 9 — a 3 sin 9 — a 4 sin 9 = 
ati sin 9 + a 2 sin 9 — a 3 cos 9 — a 4 cos 9 = 


J_dp 
p(? dx 

1 dp 
pc* dy 


1 du 

ax cos 2 9 — ar 2 cos 2 9 + a 3 sin 2 9 — a 4 sin 2 9 = 

c dx 


c*i sin 9 cos 9 — ar 2 sin 0 cos 9 — a 3 sin 9 cos 0 + ar 4 sin 9 cos 0 — 


u> _ 1 du 

2c c dy 


Ct? 1 

<*1 sin 9 cos 0 — a 2 sin 9 cos 0 — a 3 sin 0 cos 9 + a 4 sin 0 cos 0 H = 

2c cdx 

1 

ai sin 2 9 — 012 sin 2 9 + a 3 cos 2 9 — a 4 cos 2 9 — - — 

c dy 

1 . $ 

»i cos 9 + £*2 cos 9 — 0(3 sin 9 — a 4 sin 9 + /? cos <f> = — — 

p ox 


ai sin 0 + a 2 sin 9 + a 3 cos 5 + ar 4 cos 9 + /? cos = 


1 dp 
pdy 


This system can be solved for the strengths and angles of the waves, giving 
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Thus a scheme based on this model would, at each iteration, calculate the derviatives VU, 
and use the above formulas to solve for the unknowns ori, 02 , 0 * 3 , or 4 ,/?, u>, 0 , and <f>. Practical 
schemes based on this composition have yet to be developed, but work towards this end 
has been done by Roe [3], Deconinck [ 6 ] and Kroner [7]. This is an eight- wave model, 
decomposing the difference between two states into eight separate waves. The flux in a 
standard grid-aligned Roe scheme is based on a four wave model. 

In the second method, an approximate diagonalization of the Euler equations is sought. 
It is of the form 

+ Dx~ f Dy — S , (27a) 
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where W* is a vector of convected quantities (entropy, a component of velocity, and two 
acoustic-like variables), D x and D y are diagonal matrices of convection speeds, and S is a 
source term. In this approach, two /c-vectors are chosen in such a way as to minimize the 
components of the source term S. They are kT ^ for t he velocity-component convection and 
for the acoustic-like convection. With the two k vectors, 
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The right eigenvector matrix P* for the transformation 
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is given, in columns, by 


S = < 


P* (,) _ 


TT 


2 cK K > K 


2cK' 


2cK y 


, , j 

A 

2c 

- 

W ' 

1 

• K.W 

- ckW, 

w 

-(hoK^ •«(») 

- CU • K { 


„( J ) 


0 

P"j 2) 


-srfSr 
( 2 ) ( 2 ) 
PU/tj — pVfti; 

.«(*> 


(27c) 


(»d) 


(27.) 


(28) 


(29a) 


(29b) 


16 



1 


p* (3) = 


».«) 


u 


U 2 + V 2 
2 


±. 

2c 


(29c) 


2cK 


bf;~r> r ( ulc(1) - k(3) + c # 4 1) ) 

(i) P ' K (»y { vk(i) * k(3) + c ' c v ) ) 
(if-W ( hoKil) ,#c<3) + cu '* (l) ) 


2cK 


(29d) 


L 

These equations are general, holding for any choices of and k^K Hirsch et al [4] 
show that, in order to minimize the source terms, one needs a that is aligned locally 
with the pressure gradient, and a that is related to the strain-rate tensor. That is, k W 
is given by 

s <0 = Vp 
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and is computed from the velocity derivatives in the following way: if 
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then the propagation angle 
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is calculated; otherwise, the possible propagation angles are given by 
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The value of is then 

= cos 0e x + sin 9e y . (30e) 

The proper branch for Equation 30d is the one that maximizes the inner product k ^ 

This decomposition, as the standard grid-aligned Roe scheme, decomposes the residual 
using a four- wave model. 
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Figure 4: 24 x 8 Grid 


2.4 Test Case for Multi-Dimensional Schemes 

Three different schemes based on the above two-dimensional concepts are outlined in the 
following sections. The first is a scheme based on a downwind distribution. This scheme 
follows from ideas of Roe [3], and has been derived and implemented by Powell and Van 
Leer [8]. Schemes based on similar ideas have been developed by Struijs and Deconinck 
[9], Hirsch and Lacor [10] and Giles et al [11]. The second is a scheme based on a rotated 
Riemann solver. This idea, which is an extension of Jameson’s work for potential flow [12], 
was first put forward by Davis [13], and has been further studied by Levy et al [14]. The 
third scheme is based on a generalized Riemann solver. The results of the three schemes 
are compared to those of a standard “grid-aligned” scheme based on Roe’s approximate 
Riemann solver. The test-case is a reflected shock problem; an incoming Mach number of 
2.9 with a turning of 10° in a 3 x 1 channel. The case was run on two grids: a 24 x 8 grid and 
a 96 x 32 grid. The grids on which the cases were run, and the baseline results for first- and 
second-order grid-aligned schemes are shown in Figures 4-13. All cases were run without 
flux-limiting. In Sections 3-5, the three methods tested are discussed in detail. 


3 A Scheme Based on Downwind Distribution 


3.1 A Scheme for the 2D Convection Equation 

The heart of the first method, a downwind distribution scheme for the Euler equations, is a 
cell- vertex scheme for the two-dimensional convection equation 


du 

dt 


du 
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— + c x— + Cj - 0 . 
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Figure 5: 96 x 32 Grid 
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Figure 6: First-order Grid-Aligned Scheme — 24 x 8 Grid 
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Moo 5= 2.90 10 degrees turning 



Figure 7: First-order Grid- Aligned Scheme — 24 x 8 Grid 
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Figure 8: First-order Grid- Aligned Scheme — 96 x 32 Grid 
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Figure 9: First-order Grid- Aligned Scheme — 96 x 32 Grid 
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Figure 10: Second-order Grid-Aligned Scheme — 24 x 8 Grid 
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Figure 11: Second-order Grid- Aligned Scheme — 24 x 8 Grid 
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Figure 12: Second-order Grid- Aligned Scheme — 96 x 32 Grid 






Figure 13: Second-order Grid- Aligned Scheme — 96 x 32 Grid 


This scheme is analyzed below. On a uniform Cartesian grid, the residual for the convection 
equation is given by 
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These cell-face averages may be approximated using trapezoidal integration, giving the 
second-o r der appioximation 


where the semi-integer index denotes a average over a cell-face, i.e. 
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Using these formulas for the cell-face averages, the Fourier footprint of the flux integration 
for a cell is 


T (At Res) = —2 i i' x sin —■ cos ~ -f u y sin ^ cos ^ , 
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where the /?’s are the Fourier variables and the i/’s are the Courant numbers 
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To update the nodes, the cell-centered residual, given in Equation 32b, multiplied by At, 
will be sent to the nodes (i,j), ( i + 1 ,j), (i + 1 ,j + 1) and (i,j + 1) with weights u aw , u> ae , 
u ) ne , and u nw respectively (see Figure 14). The Fourier footprint of this distribution step is 
given by 
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, x . fix + Ay 
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(33d) 


If a simple forward-Euler time-stepping scheme is used, the net amplification factor for 
the entire scheme is 

G(v x , v y , /? r , Ay) = 1 + T (At Res) T (Dist) (33e) 

The appropriate values for the u>'s remain to be determined from the stability analysis. 
Since they correspond to convection directions, they should be determined by enforcing 
stability for the long waves (0 x ,/3 y — > 0). 

Taking the limit of Equation 33e as Ar> Ay — ♦ 0 yields the constraint 

k 2 (v x Px + VyPy) = (w 3to - w„ e ) — + (u: nw - u>«) — - — , (34a) 
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where k is some real constant. An added constraint (conservation) is that the entire residual 
must be distributed, 

u> ne + u> sw + u nw + u.’ ae = 1 . (34b) 

Also, by symmetry, if the u x and u y are such that convection is directly towards one node, 
all of the residual is distributed to that node, i.e. 


u ne = \ if u x = u y > 0 ; 
u sw = 1 if v x = u y < 0 ; 
0J nw = 1 if ~v x = u y > 0 ; 


u, e = 1 if -v x = Uy < 0 


Combining these conditions gives the distribution coefficients 

max (d, 


= 


= max 


u nw = max 


u) se = max 
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(34c) 

(34d) 

(34e) 

(34f) 


(35a) 

(35b) 

(35c) 

(35d) 


These formulas state that the residual is sent only to the nodes that define the downwind 
face, and is distributed in a weighted manner between the two nodes on that face. For a plane 
wave moving in one of the coordinate directions, the two downwind weights are equal, and 
the scheme reduces to the standard one-dimensional upwind scheme. The stability constraint 
of the scheme is 

u 2 x + v\ < 1 . (36) 

The locus of the scheme (i.e. the Fourier footprint of F(AtRes)F(Dist)) is shown for 
0°, 20° and 45° convection in Figures 15-17. The plots are generated by varying /?* and ft y 
discretely, which leads to a mesh of points within the continuous footprint of the locus. The 
circular stability boundary of forward-Euler time-stepping is circumscribed about the loci 
for reference. The loci art; very different from those of first-order upwind or central-difference 
schemes. It is the wave that is convected at 45° that, is damped the best, while waves at 0° 
(or 90°) are not damped well. This can be seen clearly in the contours of the amplification 
factor for the scheme, shown in Figures 18-20. 

Some numerical results for a convected Gaussian on a 32x32 grid are given in Figures 21- 
23. The convection directions are 0°, 20° and 45°. In each case the Gaussian propagates 
across the grid virtually undistorted. The onset of a zebra instability can be seen in the 
0° case, as predicted in the stability analysis. The amplitude of oscillations in this case is 
very small (on the order of 10 -4 ). The convergence history for each of the cases is shown in 
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Figure 15: Fourier Footprint of Scheme — 0° wave 
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Figure 16: Fourier Footprint of Scheme — 20° wave 
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Figure 17: Fourier Footprint of Scheme — 45° wave 
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Figure 18: Amplification Factor of Scheme 0° wave 
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Figure 19: Amplification Factor of Scheme — 20° wave 
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Figure 20: Amplification Factor of Scheme — 45° wave 
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Figure 21: Gaussian Convicted at 0° 


Figures 24 26. The Gaussian convects at almost one roll per iteration, so that the slope of the 
residual curve changes drastically after approximately forty iterations. The 45° case, which 
has the best high-frequency damping, converges very quickly, while the 0° case converges 
very slowly. 

The convergence characteristics of this scheme can be changed by an interesting variation 
in the way in which the downwind distribution step is implemented. In the above scheme, 
the residual for a cell was calculated, and split into contributions for each node of the cell. 
That is, a fraction cv n€ was distributed to the northeast node, a fraction uf nw to the noithwcst 
node, and so on. A variation of this is not to split tl e residual at all, but to distribute the 
entire residual to one node, using thotcTs as probabilities. That is, the residual is distributed 
to the northeast node with probability to the northwest node with probability u) nnn and 
so on. This scheme gives the same steady-state solutions as the residual-splitting scheme, 
but the convergence characteristics of the 0° and 20° schemes are substantially different 
(Figures 27 and 28). The 45° case is not affected by he “random distribution” method. 

Apparently the randomness of the distribution step helps kill high-frequency waves. How- 
ever, in the early parts of the convergence history, this method can actually increase the 
residual (Figure 27). The best scheme is probably a hybrid of the residual-split method and 
the “random distribution’' method. 

3.2 Scheme for the Euler Equations 

Just as in the scheme for the convection equation, the scheme for the Fulei equations is 
made up of two primary steps: 

1. a residual calculation based on a flux integral; 

2. a residual distribution. 
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Figure 22: Gaussian Convected at 20° 
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Figure 23: Gaussian Convected at 45° 
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Figure 26: Convergence History for 45° Case 


v x = 0.95 Vy = 0.00 i = 0,25 



Figure 27: Convergence History for 0° “Random Distribution” Case 
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Figure 28: Convergence History for 20° “Random Distribution” Case 


Each of these steps is somewhat more complicated for the system, however. 

For the Euler equations, the residual for a cell is given by 

Iles 1+ i J+i = -i £ [Fd» - Gdx] , (37) 

where A is the cell area and the integral is taken along the cell’s boundary dA . The cells 
are now quadrilaterals, with faces lying along curvilinear coordinate lines £ = and rj = rjj. 
The boundary integral of Equation 37 is composed of contributions of fluxes normal to 
cell faces. To extend the approximation of Equations 32e and 32f so that they apply to 
the above residua], it is necessary to convert these to flux approximations. Equation 32f, 
for instance, when multiplied by c x Ay } becomes an expression for the total flux across the 
cell-face centered at (?!, j + 1 /2): 


C X U{ j C X H{ a 

c x u UJ +iAy = ^—Ay . 


(38) 


With regard to the Euler equations, this translates directly into 
in which the following notation is used: 




X i,j + 1 — - r >,j 

Vi,j + 1 — Vi,} 


\y 

At/ 


Vj+i~ Vj 
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(391,) 

(39c) 

(39,1) 




The quantity F denotes the flux normal to a cell-face, scaled such that 

Fdr] = Fdy — Gdx . 

The analog of Equation 32e is 

G .+ ijAf = o y o A « x 


with 


(39e) 

(40a) 

(40b) 
(40c) 
(40d) 

Gd( = Fdy - Gdx . (40e) 

Since, in all but the simplest cases, the grid cells are not Cartesian, the Cartesian Courant 
numbers, v x and v y must be replaced with “contravariant Courant numbers.” The contravari- 
ant Courant numbers and are related to the wave speeds normal to the cell faces. They 
are given by the formulas 


\x = x i+hj - x itj 
= y.'+ij - Vi,j 

A£ = 6+1 -6 


and 


A* 


with 


( I/ *)»'+i.j+£ — a ( c x\y c K^i) i+ i J+ i , 


(A»j®) l+ i j + x — 2 ( x, .i+i X, ’J x *+i,i+i x i+i,j) 

(&r,y) i+ i i+ 1 = o (y'd+i — y*j + — Vi+id) > 


and 


with 


At 


1 — A (cjA^j/ c y A^a:) |+ jL^ + i , 




(A^x)- + ij + i — „ ( 2 ’t+i,j x i,j T ^i.+ij) 


(41a) 

(41b) 

(41c) 

(42a) 

(42b) 

(42c) 


(A^j/) f+ i j + | — - (yi+i,j - yi,j + y.+x,i+i - j/.j+i) • 

The distribution step requires, in each cell, projection of the residual onto the columns 
of the matrix P*, giving weights r^, and multiplication of each of the resulting vectors by 
an appropriate time-step: 


k= 1 '+2-J+2 

(43a) 

( A ' { ‘) r “)P' ( ‘)) j+4 , + l - 

(43b) 
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Figure 29: Downwind Distribution Scheme — 24 x 8 Grid 


where k varies from one to four. In the above, P*^ denotes the k ih column of P*, corre- 
sponding to the A? th wave, and is the portion of time change in the state vector caused 

by the & th wave. Each may be divided between the two vertices of the downwind 

cell-face, according to the weights of Equations 35a 35d, with v x and v v replaced by ^ ^ and 


v. 


■(*) 


e.g. 


= max 



(44) 


The wave speeds and arising in the calculation of and are the diagonal 
elements of D x and D yy given in Section 2. 

Results for this scheme arc shown in Figures 29- 32. The checkerboard instability permit- 
ted by the scheme can be seen in the odd-even decoupling that takes place upstream of the 
first shock. It was necessary to damp this mode with a nonlinear artificial viscosity of the 
type developed by Jameson et al [15]. To avoid random directions k ^ and in regions of 
small velocity and pressure gradients, both were taken to be along the flow direction in these 
regions. Despite this, the k vectors are quite noisy, as can be seen in Figures 33 and 34. The 
overall results for this scheme are disappointing, in that they are almost indistinguishable 
from the results for second-order grid-aligned upwinding. 


4 A Scheme Based on a Rotated Riemann Problem 


In the second method, tl c pressure gradient angle 0.,, given by 


0 p — tan 1 


dp / <)y 
dpj Ox ’ 


(45a) 
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Figure 30: Downwind Distribution Scheme — 24 x 8 Grid 
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Figure 31: Downwind Distribution Scheme — 96 x 32 Grid 
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Figure 32: Downwind Distribution Scheme — 96 x 32 Grid 
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Figure 33: Downwind Distribution Scheme — Vectors 
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Figure 34: Downwind Distribution Scheme — Vectors 


is used to construct the dominant-direction angle. To avoid random angles in areas with small 
pressure gradients, the dominant-direction angle used is a blend of the pressure gradient angle 
with the flow angle, 0,. The blending is accomplished by 


cos 0 = 


cos 6 P -|- (1 - e) ^ cos 0„ 

^ + (1-0“ 


(45b) 


where a typical value is e = 0.9. 

The generic fomula for a flux normal to a cell face is 


F = F(U l ,U r ); 


(46) 


in a conventional first-order-scheme the left and right input states simply are the average 
values in the adjacent cells. Convection speeds are based on velocity components parallel 
and perpendicular to the cell face; the dominant direction, i.e., the direction in which the 
Riemann problem is solved, is normal to the cell face. 

When a dominant direction is chosen that is not normal to the cell face, the flux compo- 
nent in this direction is 

Fx = F x (U Li ,U R x ). (47) 

This flux is computed as before, with two changes. The velocity components are now in the 
new reference frame (that is, the angle 0 g is replaced with the dominant- direction angle, 0), 
and the left and right states are interpolated from the values in the nearest cell centers, as 
shown in Figure 35. The subscript _L is chosen to indicate that the direction, although not 
normal to a cell face, may be normal to another important line, such as a shock front. For 
higher-order interpolation, a second “outer” ring of cells is also used, as shown in Figure 36. 
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Figure 37: First-Order Scheme with Rotated Riemann Solver — 24 x 8 Grid 


Next, a flux normal to this dominant direction, 

F|| “ (Ul||,Urj|) , (48) 

must be calculated; the input states for this flux are also found by interpolation. 

Note that this flux component also exists in the conventional formulation, but does not 
contribute to the flux through the cell face. 

The flux normal to the cell face is constructed by rotating the above fluxes back to the 
coordinate frame normal to the cell face: 

F (u La ,U Rx ,U L|| ,U R|| ) = cos {6 - 0 g ) Fj. - sin (0 - 0 g ) F () . (49) 

Here, 0 is the dominant-direction angle and 6 g is the angle of the normal to the cell face (the 
“grid angle”). 

Because of the necessity of both a ^parallel” and a “perpendicular” flux function, this 
scheme corresponds to using an eight-wave model for the fluid interaction at each interface. 

Results for the first-order scheme are shown in Figures 37-40. Results for the second- 
order scheme are shown in Figures 41—44. Comparison with the results of the grid-aligned 
scheme (Figures 6-13) shows a substantial improvement for the first-order scheme, but a less 
noticeable improvement for the second-order scheme. The dominant-direction vectors are 
shown in figure 45. Away from the shocks, the angle is in the flow direction. In the vicinity 
of the shocks, it is normal to the shocks. 
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Figure 38: First-Order Scheme with Rotated Riemann Solver — 24 x 8 Grid 
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Figure 39: First-Order Scheme with Rotated Riemann Solver — 96 x 32 Grid 
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Figure 40: First-Order Scheme with Rotated Riemann Solver — 96 X 32 Grid 
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Figure 41: Second-Order Scheme with Rotated Riemann Solver — 24 x 8 Grid 
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Figure 42: Second-Order Scheme with Rotated Riemann Solver — 24 x 8 Grid 
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Figure 43: Second- Order Scheme with Rotated Riemann Solver — 96 x 32 Grid 
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Figure 44: Second-Order Scheme with Rotated Riemann Solver — 96 x 32 Grid 
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Figure 45: Second-Order Scheme with Rotated Riemann Solver — Upwinding Directions 
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Figure 46: Geometric Interpretation of Generalized Riemann Solver 

5 A Scheme Based on a Generalized Riemann Solver 

The third scheme is based on the velocity- difference direction 


9 = tan 1 


f VR - V_ L 

\UR - u L 



(50) 


Depending on the relative magnitudes of the pressure difference and velocity difference, the 
difference between the two states Ul and Ur, is broken down into different combinations of 
waves: 


• if 


A p > pc\J (Au) 2 4- (Aw) 2 , 
two acoustic waves and one entropy wave all running in the ^-direction; 


• if 


A p < pc\J{Kuf 4- (An) 2 , 


one acoustic wave and one entropy wave running in the ^-direction, plus one shear 
wave running in the direction 9 4- tt/2. 


The idea behind this approach can be seen graphically in Figure 46. The axes in the 
figure are the velocity differences Au and Au, and the pressure difference Ap. The left 
state, Ul, is at the origin. For a right state Ur, inside the cone, such as Urj, the pressure 
difference is dominant, and the difference between Ul an< f Ur is described as two acoustic 
waves (the entropy wave, corresponding to a difference in density between _the two states, is 
not shown in the figure). For a right state Ur outside the cone, such as Ur^ or Ur^, the 
difference between the left and right states is described by one acoustic wave and one shear 
wave. The shear waves show up as lines in the (Au, Au) plane (no pressure difference), the 
acoustic waves show up as lines parallel to the rays generating the cone. 

In practice the scheme has to be slightly elaborated, because nonlinear feedback makes 
it necessary to ^freeze” the angle 9 over many interations. If the latest frozen value is 9\, 
an additional weak shear wave connects Ur to the plane representing the ^-direction. This 


45 



shear wave moves in the #1 direction; its graphical representation is a line segment in the 
direction 0 1 + tt/2, connecting the 9 - plane in the figure to the “frozen” 0 t plane. 

The first step of the scheme is to determine which set of waves to use. In the case 


-r| < W(Au) 2 + (An) 2 cos ( 0 - 0 X )| , 


(51a) 


one acoustic wave, two shear waves and one entropy wave are used. In the sub-case 

Apcos(0 — Ox) < 0 (51b) 

an acoustic wave of the first type is used; in the sub-case 

Apcos (0 — 9\) > 0 ( 51 c) 

an acoustic wave of the second type is used. In the case 

J§ ^ \\f ( Au ) 2 + ( Au ) 2 cos (0-0 1) , (51d) 

two acoustic waves, a shear wave and an entropy wave are used. 

The interface fluxes are then calculated as for Roe’s scheme, with 

F(U l , U r ) = | (F l + F r ) - ^ p |A fc | AV k R k , (52) 

with the definitions of the wave parameters depending on which model is used. 

For the case of two acoustic waves, a shear wave and an entropy wave, the eigenvectors 
are given by 


(51c) 


(51d) 


R,= 


U — C COS 0 ! 

v — csin 0 1 

h 0 — c(u cos 0i + v sin ) 


(53a) 


IU = 


— sin 0\ 
cos 6] 

(u cos#! — usin#!) 


(53b) 


(53c) 


\ {u 2 + v 2 ) 
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1 


R.4 = 


U + C COS 0 ( 
v + csin 0i 


(53d) 


^ h 0 + c(u cos Oi + v sin &i) 
The wave strengths for this case are given by 


A Vi = ^ P ~ pc\/(Au) 2 + I An) 2 cos ( 9 - 0\ 

AV 2 = p\J (An) 2 + (An) 2 sin ( 0 - 0i) 

AV 3 = A p- -^Ap 


AV 4 = ^ (Ap + pc ^ ( Au) 2 + (An) 2 cos (0 - 0i )) 
and the wave speeds by 

A, = (ti cos 0i + v sin 0 X — c) cos (0i — 0 ff ) 

A 2 = (ucos0i + u sin ^x) cos (0x — 0 g ) 

A;i = (u COS 01 + t)sin0i)cos(0x — &g) 

A 4 = (u cos 0i + v sin 0i + d) cos (0x — 0 g ) 

where 0 g is the “grid angle,” i.e. the angle of the normal to the interface, 
the projection onto the interface normal has been taken into account. 

For the case of an acoustic wave of the first type, two shear waves and 
the vectors are given by 


(54a) 

(54b) 

(54c) 


(54d) 


(55a) 

(55b) 

(55c) 

(55d) 

In these speeds, 
an entropy wave, 


Ri = 


tt — C COS 0 i 
v — csin 0] 

h 0 — c (it cos 0i + v sin 0 } ) 


R, = 


— sin 0i 
cos 0i 


(v cos 0i — u sin 0i) 


47 


r 3 


1 

u 

V 


( 


R4 


The wave i 


v 


V H“ 2 + ” 2 ) ) 
o 

— cos 9 \ 

— sin 9 \ 

(Osin 0j + 0 cos 9 \ ) 


\ 


: given by 


) 


For the case of an acoustic wave 
wave, the vectors are given by 


/ 


Ri 


V 


of the second type, 

0 

— cos 9 X 

— sin 9 1 

— (Osin 0 j -f 0 cos 9i) 


(56c) 


(56d) 


A Vi = 

Ap 
c 2 

(57a) 

af 2 - 

p\J (Aw ) 2 + (Aw ) 2 sin (9 — 0 X ) 

(57b) 

av 3 = 

A P - 

(57c) 

av 4 = 

t Ap — pc\j (Aw ) 2 + (Aw ) 2 cos (9 — 0 X )^ 

(57d) 

and the wave speeds by 



Ai 

= (0 cos 0 X + 0 sin 0 X — c) cos (0 X — 9 g ) 

(58a) 

^2 

= (0 cos 9i -f 0 sin 0 X ) cos (0 X — 9 g ) 

(58b) 

A 3 

= ( u cos 9\ + 0 sin 9i ) cos (0 X — 9 g ) 

(58c) 

A 4 

= (Osin — 0 cos 0i ) sin (0i — 0 g ) . 

(58d) 


and an 


/ 
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0 

— sin Oi 
cos 61 

(v cos 6 \ — usin $i) 

1 

it 

v 

i( u 2 + i> 2 ) 

1 

it + c cos Ox 
v + csin^i 

h 0 4- c (it cos Ox + t) sin 0i) 
The wave strengths for this case are given by 



and the wave speeds by 


A, = (usin0i - n cos 0i ) sin (0i — 0 g ) ( 61a ) 

A-, = (ti cos Ox + v sin 0 \ ) cos ( 0 \ — 0 g ) (616) 

A :( = (« cos Ox + t) sin ) cos (#i — 0 a ) (61c) 

A,, = (it cos Ox + v sin Ox 4- <’) cos (Ox — 0 g ) . (61 d) 


As in the standard grid-aligned scheme, the difference between two states is described in 
terms of four waves. In this scheme, however, the choice of waves depends on the relative 
magnitudes of the pressure difference and velocity difference between the states. 
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Figure 47: First-Order Scheme with Generalized Riemann Solver — 24 x 8 Grid 


Results for this scheme are shown in Figures 47-50. They show a substantial improvement 
over the first-order grid-aligned results (Figures 6-9), and a slight improvement over the 
second-order grid-aligned results (Figures 6-9). 

Further development of the generalized Riemann solver must be aimed at making it 
more robust, without giving up its resolving power. A promising approach, currently being 
investigated, is to anchor the wave decomposition to the frame of the local flow direction; 
this is accomplished by simply taking 6 \ equal to the flow angle. For normal, or almost 
normal, shocks, the results remain grid-independent; for more oblique shocks, resolution is 
lost because the flow direction is not a good indicator of the shock normal, and some weight 
should be given to the velocity-difference direction. This can be done by explaining part 
of the jump between U/, and U/j by the generalized Riemann solver based on the velocity- 
difference direction, and the remainder by the generalized Riemann solver based on the 
streamwise direction. The full flux calculation then involved seven waves. 


6 An Outlook for Multi-Dimensional Upwind Schemes 

I he outlook for genuinely multi dimensional upwind schemes is still clouded. For the re- 
flected shock cases presented here, all schemes tested showed improvements over first-order 
grid-aligned schemes, but the improvements over second-order grid-aligned were very slight. 
It is not clear whether this result is a condemnation of multi-dimensional schemes, or a con- 
firmation of grid-aligned schemes. The improvements over grid-aligned upwinding do come 
at a cost of complexity and convergence. 

The downwind distribution scheme fundamentally differs from the other two schemes 
presented here: it is a cell-vertex scheme. This implies better accuracy on stretched meshes, 
and a simpler extension to unstructured meshes than for the cell-centered schemes. However, 
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Figure 48: First-Order Scheme with Generalized Riemann Solver — 24 x 8 Grid 


Moo = 2.90 10 degrees turning 
Density Contours 

1 0.9000 

2 1.0000 

3 1.1000 

4 1.2000 

5 1.3000 

6 1.4000 

7 1.5000 

8 1.6000 

9 1.7000 

10 1.8000 

11 1.9000 

12 2.0000 

13 2.1000 

14 2.2000 

15 2.3000 

16 2.4000 

17 2.5000 

18 2.6000 

19 2.7000 

20 2.8000 


x 



Figure 49: First-Order Scheme with Generalized Riemann Solver — 96 x 32 Grid 
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Moo — 2.90 10 degrees turning 



Figure 50: First-Order Scheme with Generalized Riemann Solver — 96 x 32 Grid 


this scheme is also the least diffusive of the three, and the only one requiring additional 
damping to reach a converged steady state. There are also some difficulties associated with 
propagation directions that lie along grid lines. The random-distribution method presented 
here might alleviate these problems, but it relies on the fact that the cell-centered residuals, 
and not just the changes at the nodes, are zero when the steady state is reached. This may 
not be generally true; non-zero cell residuals could add up, in the distribution step, to give 
zero changes at the nodes. If this is the case, the random distribution scheme might not lead 
to convergence, or to the proper steady state. 

The scheme based on a rotated Riemann solver is the most robust of the three schemes. 
Its cost, however, is fairly high. At each interface, almost twice the work of a grid-aligned 
Riemann solver is necessary to calculate the fluxes. Furthermore, the interpolation of the 
state vectors has an associated cost, both in computational work and in added storage. In 
three dimensions, the added storage could be prohibitive. The interpolation technique would 
be cumbersome on an unstructured grid. 

The scheme based on a generalized Riemann solver seems, at this point in time, the most 
promising of the three approaches. The computational work for the scheme as outlined is 
only slightly more than that of a grid-aligned scheme. The nonlinearity added by basing 
the scheme on an angle calculated from local flow variables is detrimental to convergence, 
but ways to get around this are being developed. In one approach, a portion of the differ- 
ence between the left and right states (approximately one-half) is decomposed according to 
the generalized Riemann solver, and the remainder is decomposed according to a stream- 
aligned Riemann solver. This approach, however, adds to the cost of the flux computation, 
approximately doubling it. 

Clearly, the simple shock-reflection problem used in this paper cannot truly test these 
multi-dimensional schemes. More sophisticated problems, particularly ones in which shear 
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layers lie oblique to the grid, must be studied. Also, the possible range of schemes based on 
multi- dimensional ideas has not nearly been covered yet. For a topic in which research has 
begun only recently, though, the results are promising. 
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